Thermomechanical properties of a single hexagonal boron nitride sheet
m
o
(N
Oh!
<
(N
(N
a
o
>
^'
o
m
X
Sandeep Kumar Singh^, M. Neek-Amal^^^*, S. Costamagna^"'^, F. M. Peeters ^
^Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerpen, Belgium
'^Department of Physics, Shahid Rajaee University, Lavizan, Tehran 16785-136, Iran
^Facultad de Ciencias Exactas Ingenieria y Agrimensura,
Universidad Nacional de Rosario and Instituto de Fisica Rosario,
Bv. 27 de Fehrero 210 his, 2000 Rosario, Argentina.
(Dated: April 23, 2013)
Using atomistic simulations we investigate the thermodynamical properties of a single atomic layer
of hexagonal boron nitride (h-BN). The thermal induced ripples, heat capacity, and thermal lattice
expansion of large scale h-BN sheets are determined and compared to those found for graphene
(GE) for temperatures up to 1000 K. By analyzing the mean square height fluctuations {h'^) and
the height-height correlation function H{q) we found that the h-BN sheet is a less stiff material as
compared to graphene. The bending rigidity of h-BN: i) is about 16% smaller than the one of GE at
room temperature (300 K), and ii) increases with temperature as in GE. The difference in stiffness
between h-BN and GE results in unequal responses to external uniaxial and shear stress and different
buckling transitions. In contrast to a GE sheet, the buckling transition of a h-BN sheet depends
strongly on the direction of the applied compression. The molar heat capacity, thermal expansion
coefficient and the Gruneisen parameter are estimated to be 25.2 JmoP^K"'^, 7.2xlO~®K~^ and
0.89, respectively.
I. INTRODUCTION
A single layer of hexagonal boron-nitride (h-BN) is
a wide gap insulator that is very promising for opto-
electronic technologies [l|, 0], tunnel devices and field-
effect transistors [SHa]. According to the well known
Mermin- Wagner theorem Q, the stability of any two
dimensional crystal is only possible in the presence of
atomic corrugations which distort the perfect honey-
comb lattice and allow the atoms to explore the out-of-
plane direction. Experimental observations have found
ripples in suspended sheets of GE [3, Q and atomistic
simulations suggest that the strong bonds between the
carbon atoms determine the features of these ripples [9j .
Understanding the behavior of the ripples is important
for many reasons [Ifll- They affect the electronic trans-
port properties, e. g., in GE the ripples are believed to
be one of the dominant scattering sources which limits
the electron mobility [ll|, [l2] ■
h-BN ribbons doped by carbon has recently been pro-
posed [3, [l^. In addition, BN based nanostructures
are potential materials for thermal management applica-
tions |14l4l3 | because of their high thermal conductivity
and sensitivity to isotopic substitution and etc. There-
fore, the knowledge of the shape and the temperature
dependence of the intrinsic ripples is fundamental to de-
vise novel nano-devices based on this material.
Both GE and h-BN sheets have a honeycomb lat-
tice structure, however the different electronic proper-
ties and phonon band structure |20l - |22| results in un-
equal morphologies and corrugations. Transmission elec-
tron microscopy is widely used to resolve the individ-
* corresponding author: ncckanial@srttu.edu
Armchair
O Boron O Nitrogen
FIG. 1: (Color online) Schematic view of the single h-BN
sheet. Smaller-yellow (bigger-blue) circles refer to the B (N)
atoms. The rectangles indicate the atoms that are fixed dur-
ing compression. Dashed (straight) lines correspond to arm-
chair (zig-zag) uniaxial compression in the direction given by
the arrows. Open arrows indicate the shear stress applied in
the armchair direction.
ual atoms in suspended h-BN sheets [23| where ripples
inherently exist. First-principle calculations have been
performed using small unit cells, periodically replicated,
which are unable to model long wavelength corrugations
which require thousands of atoms [2J| while the mechan-
ical properties of a h-BN sheet can be estimated by us-
ing DFT 123]. Conversely, atomistic simulations using
molecular dynamics simulations (MD) enable to study
thermo-mechanical properties directly on large scale sam-
ples. The modified TersofF potential [26,] (TP) (parame-
terized originally for carbon and silicon) with various set
of parameters have shown to predict reasonable values for
the thermo- mechanical properties of the h-BN sheet. In
the pioneer work by Albe et al re-parameterized TP was
used to study the impact of energetic boron and nitrogen
atoms on a cubic-BN target 27]. Some other groups have
also parameterized TP using various experimental data,
e.g. Sekkal et al [2a] treated h-BN as a one-component
system, using the same potential parameters for both
boron (B) and nitrogen (N) (neglecting the B-N interac-
tion) to investigate the structural properties. Matsunaga
et al Q proposed the TP of B in order to simulate cu-
bic boron carbonitrides which are cornpatible with the
parameters of nitrogen fitted by Kroll [30| , and recently,
Liao et al [31| and Sevik et al [32| reported TP param-
eters that were used to study the mechanical properties
and the thermal conductivity of a h-BN sheet, respec-
tively.
In this study, we investigate the thermal rippling be-
havior of free standing monolayer h-BN by using state
of the art molecular dynamics (MD) simulations. We
adopted the TP potential re-parametrized by Sevik et
al which has been shown to represent the experimen-
tal structure and the phonon dispersion of the two-
dimensional h-BN sheet. We found that h-BN is more
corrugated and a less stiff material as compared to GE.
The height-height correlations can be explained by the
theory for continuum membranes [33]. In addition, we
report results of both uniaxial and shear stress of a h-BN
sheet and compare it with those found for GE. The buck-
ling transition for compressed h-BN occurs earlier than
for GE and the induced pattern of ripples when subjected
to either uniaxial or shear stress shows significant differ-
ences.
This paper is organized as follows. In Sec. II, we in-
troduce the atomistic model and the simulation method.
Then, in Sec. Ill we analyze the behavior of the thermal
ripples of a h-BN sheet. Here, we obtain the bending
rigidity, the heat capacity, the thermal expansion coeffi-
cients, and we study the buckling transition of the h-BN
sheet when uniaxial external strain and shear stress are
applied. All the results are compared with the ones ob-
tained for a GE sheet. Finally, in Sec. IV, we present
our conclusions.
II. METHODS
Classical atomistic molecular dynamics (MD) simula-
tion is employed to simulate large scale samples of h-BN
and GE at temperatures varying from 10 K up to 1000 K.
We used a modified Tersoff potential (which is an ordi-
nary defined potential in the LAMMPS package [33.l35|)
according to the parameters proposed by Sevik et al for
the h-BN sheet. All the parameters and a detailed de-
scription of the potential energy is given in Ref . [32| . To
simulate GE we have used the AIREBO potential [3g
which is suitable for simulating hydrocarbons. We em-
ployed NPT ensemble simulations with P=0 using the
Nose-Hoover thermostat which enables us to allow the
size of the system to equilibrate (i.e. the size of the
system is not fixed). All the reported values have been
computed averaging over 300 — 400 snapshots taken over
uncorrelated samples.
We start with a square shaped h-BN sheet with
periodic boundary conditions and initial dimensions
322A x32ll (315A xSlsA for GE) in the x and y direc-
tion, respectively, which correspond to a total number
of A^=37888 atoms, and which are sufficiently large in
order to capture the long-wavelength regime. Periodic
boundary conditions were adopted in both directions.
To analyze the thermal ripples we computed the
diffraction pattern which is typically studied in experi-
ments to detect the size and shape of the corrugations [3| .
We obtained also the mean square value of the out-of
plane displacements JJh^) of the atoms and, by follow-
ing previous works [33, |3g], the height- height correla-
tion function {H(q)) which was determined by consid-
ering an average of the height between the first neigh-
bors of each atom. The latter was shown to follow a
q^^ behavior that is expected from the theory of con-
tinuum two-dimensional membranes at large values of q
in the harmonic approximation (see below). To analyze
the differences between the strain induced corrugations
in the h-BN and the GE sheets we applied uniaxial and
shear stress on both systems as is schematically shown
in Fig. [T] In order to apply strain we set the force on
the two ends equal to zero and move the end atoms with
an infinitesimal compression step {dx = O.OlA) in the
desired direction. After each compression step we wait
for 2 ps to allow the system to relax and to ensure that
the system is in thermal equilibrium [39|. Uniaxial com-
pressive stress is applied along the zig-zag or armchair
direction separately, and the shear stress is applied on
the armchair edges with the opposite velocity for the top
and bottom edges. The details of the used method of ap-
plying the boundary stress can be found in our previous
studies [H-Iill
The TP function [H] used in the LAMMPS pack-
age [SJ, I33 can be expressed as
E
e^.4e
^^j)
with
i) = Yl Yl Mrrj)[fRinj) + bijfAirrj)]
(1)
(2)
j(>«)
where /c, /r and /a are cutoff functions, the repulsive
pair term, and the attractive pair term, respectively. Ty
and bij are respectively the distance from atom i to atom
j and the bond order function. The use of TP disregard
contributions coming from charge re-distribution which
may occur in an ionic crystal. The inclusion of this effect
in h-BN modifies the phonon spectrum significantly only
for energies corresponding to the optical modes [42|, |43| .
The large scale thermal ripples addressed here are gov-
erned mainly by the transversal acoustic mode (TA),
(c)
0.1
n 1 1 i —
-3 -2 -1
0.05
-0.05
-0.1
-0.1 -0.05 0.05 0.1
FIG. 2: (Color online) (a) Contour plot of the heights for an arbitrary snapshot taken during the MD simulation at 300 K.
(b) Corresponding simulated diffraction pattern for the sample shown in (a), (c) Zoom view of (b) around q=(0,0).
,\\ ' I ' I ' I
(a)
GE
h-BN
200 400 600 800 1000
T [K]
0.125 0.25 0.5
q [A"']
FIG. 3: (Color onhne) Variation of (a) < h^ > and (b)
the bending rigidity versus temperature for both h-BN (open
squares) and GE (solid circles), (c) Height-height correlation
function H{q) of h-BN at two different temperatures as is
indicated. The dashed lines are the harmonic q~* law.
which accounts for out-of plane fluctuations, and it cou-
ples with the in-plane modes which renormalizes the long
wave-length behavior, e.g. the bendiiig rigidity k can be
calculated directly from the TA mode l4J] . Therefore, the
charge redistribution is expected not to affect the ther-
mal fluctuations analyzed here and the use of the TP is
justified ^^.
III. RESULTS AND DISCUSSION
A. Thermal excited ripples
In Fig.[31Ja) we show a height contour plot of the atoms
of the h-BN sheet for an arbitrary snapshot taken during
the MD simulation at 300 K. The corresponding modeled
diffraction pattern is shown in Fig.[2][b). This pattern is
very similar to the one obtained for the GE sheet Q with
the main difference in the distance between the Bragg
points due to the unequal lattice constant of h-BN and
GE. From the zoom plot around q=(0,0) (Fig.[2][c)) one
observes the local structure of the central Bragg point
for these intrinsic thermal ripples of h-BN. Notice that
the lack of the presence of the q=(0,0) component is a
consequence of the absence of a perfectly flat h-BN sheet.
The signatures of the thermal induced ripples can also
be seen in the mean square value of the out-of-plane
fluctuations {h"^)- In Fig. ^a.) we show (/i^) as func-
tion of temperature. In comparison with GE (included
also in this figure for comparative purposes) we observe
that (/i^) is always larger for the h-BN sheet. The
weaker (stronger) the atomic B-N (C-C) bonds, the larger
(smaller) the corrugations in h-BN (GE). Notice that the
B-N bond is not pure covalent and it has an ionic charac-
ter which is due to the different electronegativity between
the two elements, i.e. the electrons are localized closer
to the N atoms rather than the B atoms. Although this
partially ionic structure of the h-BN layer increases the
interlayer interaction resulting in a larger hardness of 3D
bulk h-BN relative to graphite, it makes the single layer
of BN less stiff than graphene. Moreover, the stacking
of h-BN sheets in bulk h-BN is AAA stacking which is
different for the ABC stacking in graphite j23i] .
B. Bending rigidity K
Accordingly to the two-dimensional theory of contin-
uum membranes the height-height correlation function,
in the harmonic approximation, where the out-of-plane
and in-plane modes are decoupled, is expected to behave
as
iJ(g) = (|%)|2} =
NkuT
kSoq
A '
(3)
where k is the bending rigidity of the membrane, iV is
the number of atoms of the sample, 5*0 the surface area
per atom and fcs is the Boltzmann constant. In the large
wavelength limit, i.e. for g — > 0, the height fluctuations
are suppressed by anharmonic couplings between bend-
ing and stretching modes giving rise to a renormalized q-
dependent bending rigidity and hence Eq. ([3]) is no longer
valid [Hill 0- In Fig. ^h) we show k for h-BN cal-
culated from the harmonic part of H(q) between q=0.5
A~^ and q=l A~^. In agreement with the larger values
of (/i^), we observe that the h-BN membrane posses a
smaller k as compared to GE and in the whole temper-
ature range it is about 16 % smaller at room tempera-
ture (300 K). Note that k for both h-BN and GE increase
with temperature. In Fig. |3Kc) we show H{q) at 200 K
and 1000 K were the harmonic behavior can be clearly
observed (fitted with a dashed line) and, as expected,
with increasing temperature H{q) is shifted to larger q.
This figure also reveals that the ripples are not charac-
terized by an unique wave-length and instead follows the
behavior expected from continuum membrane theory.
Before ending this section it is worthwhile to
investigate an alternative method to estimate the
bending rigidity (flexural rigidity). A common
method for calculating the bending rigidity of
BN-sheet is by performing simulations of BN-
nanotubes as a function of its radius (R) of the
curved tubes and then extrapolating the results
to R^> oo. Hence, one can calculate the elastic en-
ergy of the nanotube as a function of the inverse
square of the radius, E ~ ^kR^^. This method was
used in Ref. [25] and in our previous work [41] to
calculate the zero temperature bending rigidity of
GE and h-BN which were found to be 1.46 eV and
1.29 eV, respectively. The result obtained with
the Tersoff potential using Eq. (3) is less than the
result of Ref. [25]. The bending rigidity of GE
is larger than h-BN with about 0.15 eV in agree-
ment with Ref. [25]. In order to have an indepen-
dent check we estimated the bending rigidity of
a BN sheet by performing several ground state
simulations for (n,n) BN-nanotubes with n=5,..,
20 using the Tersoff potential. Figure 4 shows the
variation of the strain energy per atom as func-
tion of the inverse square radius of the tube. Fit-
ting a line to the data and dividing by the area of
half of unit cell atom results in k = 0.86 eV. The
result of Tersoff potential agrees well with our
0.00
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14
FIG. 4: (Color online) Variation of strain energy versus
inverse square of BN-nanotube radius using Tersoff
potential.
10 200 400 600 800 1000
T[K]
FIG. 5: (Color online) Variation of (a) the total energy per
atom and (b) the averaged bond length versus temperature
for h-BN and GE (The error bar is less than 10~* eV/atom
for the total energy).
estimation for k using Eq. (3) (Ref. [25]). Thus
we can conclude that result based on the Tersoff
potential underestimate the bending rigidity as
compared to the result from DFT.
C. Heat capacity, thermal expansion and
Gruneisen parameter
The variation of the total energy per atom with tem-
perature (> lOK) is shown in Fig. [5ja). The total en-
ergy varies linearly with temperature and gives the corre-
sponding lattice contribution to the molar heat capacity
at constant pressure (the average size of the system after
relaxation is taken constant) which for the h-BN sheet
resuhs into
Cp =
dEr
~dT
= 25.2(3) J mor^iC-^
(4)
which is comparable to the one for a GE sheet, i.e.
24.5(9) J mol~^K~^ and close to the well known classical
molar heat capacity, i.e., C = 3K ~ 24.94 J mol~^K^^,
i.e., the Dulong-Petit value, where SR is the universal gas
constant. Notice that the heat capacity is a tem-
perature dependent parameter that will decrease
for temperature below the Debye temperature.
However, the classical MD simulation gives the
correct high temperature asymptotic limit i.e.,
the Dulong-Petit value, but fails in the low and
intermediate temperature range.
In Fig. [5jb) we show the variation of the averaged B-
N bond length with temperature. The linear behavior
enables us to calculate the linear thermal expansion of
the lattice parameter of a h-BN sheet as
1 da
a dT
lBN = -^ = 7.2(1) xlO-^K-'
(5)
where a= 1.442 A is the zero temperature lattice parame-
ter. The obtained 7 is 33% larger than the one measured
for cubic BN, i.e. 4.8xlO~^K~^ [43| and is comparable
to the one found for GE, i.e. jge = 10.0(7) x IQ-^ K-\
The latter (i.e. "(nE) is in good agreement with our pre-
vious studies |39l.l40|.
From the obtained values of 7 and C we can compute
the Gruneisen parameter
asN = -pr ^ 0.89,
Cp
(6)
where B is the two dimensional bulk modulus and p
is the mass density. Using Bh^BN=SeVk^'^ [Bqe =
12.7 eVA-2 HI)), ph-BN = 2A.82/Sh-BN - 3.81 x
lO^-^gm-^ [pQj^ = 7.6 X 10-4gm-2) for h-BN (GE), we
obtain ubn = 0.89 (aGB=1.2). Note that the estimated
value of aoE is better than the one found in Ref. |49|, i.e.
-0.2, and is closer to the experimental result, i.e. 2.0 |50| .
D. Buckling transition
The different stiffness between the h-BN and GE sheets
results in different buckling transitions. To study this, we
compressed the systems uniaxially where we considered
both zig-zag and armchair directions. This was done by
fixing one row of atoms in each edge during the compres-
sion steps as is indicated schematically in Fig. [TJ The
compression rate was taken p = 0.5m/ s which is simi-
lar to the one used in our previous study [40|, and small
enough to guarantee that the system is in equilibrium
during the whole compression process |5l| . The simula-
tions were carried out at room temperature. Figure [BJ^a)
shows the variation oi < h^ > versus strain, i.e. e = pt/l
0.01 0.1
strain (%)
(c)
150
100
50
<
-50
llHWUn r'l i 1
-100
-150
;,<
,t f ,t t ,f
1
7
6
5
4
3[A]
2
1
-1
-2
-3
-4
-100
0.
x(A)
100
FIG. 6: (Color online) Variation of < ft^ > versus external
uniaxial strain in the (a) armchair and (b) zig-zag direction.
Contour plot of the heights of compressed h-BN sheet samples
subjected to a compressive strain for fixed < h^ >= 20 A^
for compression in (c) zig-zag direction of the h-BN sheet.
Arrows indicate the stress direction.
where t is the time (after starting the compression) and
I is the initial box length in the direction of the compres-
sion. The buckling transition for the h-BN sheet kept
at room temperature is found to be 0.6 (0.1) which is
smaller than the one for GE, i.e. 1.0 (1.2), for uniaxial
compression along the zig-zag (armchair) direction 52].
Hence, when the compression is applied in the zig-zag
direction, the h-BN sheet resists much more against the
external uniaxial stress as compared to the case the stress
is applied along the armchair direction. The smaller
critical strain at which the buckling transition in
the armchair direction of the h-BN sheet takes
place indicates its more sensitive nature to exter-
nal uniaxial stress [53]. Although the same argument
holds for GE, the difference between the two directions
is much smaller. Notice that DFT calculations re-
sult in a deviation in bending rigidity (flexural
rigidity) between ZZ BN-nanotube and AC BN-
0.65
0.0001 , ,0.0002
U' (A)"'
0.0003
FIG. 7: (Color online) Variation of buckling strain with the
longitudinal length of the BN sheet which is compressed along
the zig-zag direction.
nanotube (ZZ becomes larger than AC) for radius
smaller than ~ SA [25] while for larger radius they
are the same.
A contour plot of the buckled h-BN sheet with (ft,^)=20
A^ and compressed in the zig-zag direction is shown in
Fig. IHc). The obtained buckling transition for GE, i.e.
0.8%, agrees very well with our previous studies [39| and
is in the range of the experimental measured buckling
transition for suspended GE, i.e. 0.5%-0.7% ^E^.
It is also interesting to compare the buckling strain
with those predicted by the theory of loaded beam. Euler
theory for a two end-loaded beam having length L pre-
dicts that efc oc L~^ [40, lla^. The first demonstration
by MD of Euler buckling in nanostructures goes
back to the pioneer work by Yakobson et al |.57J .
We performed several simulations in order to find the
length dependence of the buckling strain for BN sheets
which are e.g. compressed along the zig-zag direction.
Figure [7] shows the variation of ei, with the longitudinal
length (L) of the zig-zag BN sheet. The solid line is a fit
to L~^ and the symbols are the results of our MD sim-
ulations which are qualitatively in agreement with the
theory of loaded beam.
Finally, we report the results for h-BN and GE sheets
under shear stress. Here we applied the stress only in the
armchair direction as is described schematically in Fig.[T]
We found that, due to the larger stiffness, to reach the
same value of {h"^) a larger shear stress has to be applied
in GE as compared to h-BN. In Figs. |8l^a,b) one can ob-
serve significant differences between the corrugated con-
figurations of h-BN and GE sheets. These samples were
subjected to a shear stress of e = 1.5%. While h-BN is
highly sensitive to shear stress and deforms easily, GE
can resist larger stress values due to its larger rigidity.
IV. CONCLUSIONS
The thermal properties of a boron nitride sheet were
studied using large scale atomistic simulations. We
showed that the scaling properties of a h-BN sheet fol-
lows closely the results of membrane theory and hence
the thermal excited ripples are not characterized by any
particular wave-length. Using the harmonic part of the
height-height correlation function we found an increas-
ing bending rigidity with temperature which is smaller
than the one of graphene. We found that the buckling
transition for h-BN depends on the applied compression
direction and is much smaller than the one of graphene.
The obtained molar heat capacity agrees very well with
the well-known Dulong-Petit number, 25.2 J mol~^K^^
and the thermal expansion coefficient was found to be
positive and equal to 7.2 xlQ~^K~^. The Gruneisen
parameter 0.89 is found to be smaller than the one for
graphene, i.e. 1.2. We showed that the different stiffness
between the GE and h-BN sheets leads to different pat-
terns of deformations in the presence of either uniaxial
or shear stress.
Acknowledgments
We thank K. H. Michel and D. A. Kirilenko for
their useful comments on the manuscript. M. N.-A.
is supported by the EU-Marie Curie IIF postdoc Fel-
lowship/299855. S. Costamagna is supported by the
Belgian Science Foundation (BELSPO). This work is
supported by the ESF-EuroGRAPHENE project GON-
GRAN, the Flemish Science Foundation (FWO-Vl), and
the Methusalem programme of the Flemish Government.
[1] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Phys.
Rev. B 51, 6868 (1995).
[2] K. Watanabe, T. Taniguchi, and H. Kanda, Nat. Mater.
3, 404 (2004).
[3] L. BritneU, R. V. Gorbachev, R. Jalil, B. D. Belle, F.
Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A.
S. Mayorov, N. M. R. Peres, A. H. Castro Neto, J. Leist,
A. K. Geim, L. A. Ponomarenko, and K. S. Novoselov,
Nano Lett. 12, 1707 (2012).
[4] G-H. Lee, Y-J. Yu, C. Lee, C. Dean, K. L. Shepard, P.
Kim, and J. Hone, Appl. Phys. Lett. 99, 243114 (2011).
[5] J. Beheshtian, A. Sadeghi, M. Neek-Amal, K. Michel,
and F. M. Peeters, Phys. Rev. B 86, 195433 (2012).
[6] N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133
(1966).
[7] J. C. Meyer, A. Geim, M. I. Katsnelson, K. S. Novoselov,
0.8
0.6
0.4
0.2 .
[A]
-0.2
-0.4
-0.6
-0.8
FIG. 8: Effect of a shear stress of about 1.5 % on the out-of-plane deformation of (a) the h-BN sheet and (b) the graphene
sheet. Arrows indicate the stress direction on the armchair edges. Note the larger values of the heights in h-BN.
[9
[lo;
[11
[12;
[is;
[14
[is;
[16;
[it;
[is;
[19
[2o;
[21
[22;
[23;
[24
T. J. Booth, and S. Roth, Nature (London) 446 60
(2007).
D. A. Kirilenko, A. T. Dideykin, and G. Van Tendeloo,
Phys. Rev. B 84, 235417 (2011).
A. Fasohno, J. H. Los, and M. I. Katsnelson, Nat. Mater.
6, 85 (2007).
S. Costamagna, M. Neek-Amal, J. H. Los, and F. M.
Peelers, Phys. Rev. B 86, 041408(R) (2012).
M. L Katnelson, and A. K. Geim, Phil. Trans. R. Soc. A
366, 195 (2008).
S. Costamagna, O. Hernandez, and A. Dobry, Phys. Rev.
B 81, 115421 (2010).
J. Jung, Z. Qiao, Q. Niu, and A. H. MacDonald, Nano
Letters, 12, 2936 (2012).
T. Ouyang, Y. Chen, Y. Xie, K. Yang, Z. Bao, and J.
Zhong, Nanotechnology 21, 245701 (2010).
I. Savic, D. A. Stewart, and N. Mingo, Phys. Rev. B 78,
235434 (2008).
D. A. Stewart, L Savic, and N. Mingo, Nano Lett. 9, 81
(2009).
C. W. Chang, A. M. Fennimore, A. Afanasiev, D. Okawa,
T. Ikuno, H. Garcia, D. Li, A. Majumdar, and A. Zettl,
Phys. Rev. Lett. 97, 085901 (2006).
H. Shen, Comp. Mater. Sci. 47, 220 (2009).
Y. Xiao, X. H. Yan, J. X. Cao, J. W. Ding, Y. L. Mao,
and J. Xiang, Phys. Rev. B 69, 205415 (2004).
M. Sakurai, Y. Sakai, and S. Saito, J. Phys.: Conf. Ser.
302, 012018 (2011).
M. Topsakal, E. Aktiirk, and S. Ciraci, Phys. Rev. B 79,
115442 (2009).
D. Golberg, Y. Bando, Y. Huang, T. Terao, M. Mitome,
C. Tang, and C. Zhi, ACS Nano 4, 2979 (2010).
N. Alem, R. Erni, C. Kisielowski, M. D. Rossell, W. Gan-
nett, and A. Zettl, Phys. Rev. B 80, 155425 (2009).
H. Sahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E.
Aktiirk, R. T. Senger, and S. Ciraci, Phys. Rev. B 80,
155453 (2009).
[25] K. N. Kudin and G. E. Scuseria, and B. L Yakobson,
Phys. Rev. B 64, 235406 (2001).
[26] J. Tersoff, Phys. Rev. B 37, 6991 (1988).
[27] K. Albe and W. MoUer, Comput. Mater. Sci. 10, 111
(1998).
[28] W. Sekkal, B. Bouhafs, H. Aourag, and M. Cartier, J.
Phys.: Condens. Matter 10, 4975 (1998).
[29] K. Matsunaga, C. Fisher, and H. Matsubara, Japan J.
Appl. Phys. 39, 48 (2000).
[30] P. M. KroU, PhD Thesis (Technische Hochschule Darm-
stadt, 1996).
[31] M. L. Liao, Y. C. Wang, S. P. Ju, T. W. Lien, and L. F.
Huang, J. Appl. Phys. 110, 054310 (2011).
[32] C. Sevik, A. Kinaci, J. B. Haskins, and T. Qagin, Phys.
Rev. B 84, 085409 (2011).
[33] D. Nelson, T. Piran and S. Weinberg, Statistical Mechan-
ics of Membranes and Surfaces (Word Scientific, Singa-
pore, 2004).
[34
[35^
[36
[37
[38;
[39
[40
[41
[42;
|http://lammps. sandia.gov
S. J. Plimpton, J. Comp. Phys. 117, 1 (1995).
S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem.
Phys. 112, 6472 (2000).
S. Costamagna and A. Dobry, Phys. Rev. B 83, 233401
(2011).
J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Za-
kharchenko, and A. Fasolino, Phys. Rev. B 80, 121405
(2009).
M. Neek-Amal and F. M. Peelers, Phys. Rev. B 82,
085432 (2010); Appl. Phys. Lett. 97, 153118 (2010).
M. Neek-Amal and F. M. Peelers, J. Phys.: Condens.
Matter 23, 045002 (2011).
A. Lajevardipour, M. Neek-Amal, and F. M. Peelers, J.
Phys.: Condens. Matter 24, 175303 (2012).
K. H. Michel and B. Verbeck, Phys. Status Solidi B 246,
2802 (2009).
[43] K. H. Michel and B. Verbeck, Phys. Rev. B 83, 115328
(2011).
[44] L.J. Karssemeijer and A. Fasolino, Surf. Sci. 605, 1611
(2011).
[45] S. S. Han, J. K. Kang, H. M. Lee, A. C. T. van Duin, and
W. A. Goddard, J. of Chem. Phys. 123, 114703 (2005).
[46] P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett. 69,
1209 (1992).
[47] G. Demazeau, Diamond and Related Materials, 2, 197
(1993).
[48] K. V. Zakharchenko, M. L Katsnelson, and A. Fasolino,
Phys. Rev. Lett. 102, 046808 (2009).
[49] V. Perebeinos and J. Tersoff, Phys. Rev. B 79, 241409(R)
(2009).
[50] T. M. G. Mohiuddin, A. Lombardo, R. R. Nair, A.
Bonetti, G. Savini, R. Jalil, N. Bonini, D. M. Basko, C.
Galiotis, N. Marzari, K. S. Novoselov, A. K. Geim, and
A. C. Ferrari, Phys. Rev. B 79, 205433 (2009).
[51] M. Neek-Amal, Phys. Rev. Lett 106, 209701 (2011).
[52] A movie of the buckling process in a zig-zag BN sheet is
provided in the Supplementary Material.
[53] Q. Peng, W. Ji, and S. De, Comp. Mat. Science 56, 11
(2012).
[54] O. Frank, G. Tsoukleri, J. Parthenios, K. Papagelis, I.
Riaz, R. Jalil, K. S. Novoselov, and C. Galiotis, ACS
Nano 4, 3131 (2010).
[55] G. Tsoukleri, J. Parthenios, K. Papagelis, R. Jalil, A. C.
Ferrari, A. K. Geim , K. S. Novoselov, and C. Galiotis,
Small 5, 2397 (2009).
[56] J. Singer, J. Arbocz and T. Weller, Experimental Methods
in Buckling of Thin- Walled Structures (Wiley, New York,
1998). R. M. Jones, Buckling of Bars, Plates, and Shells
(Taylor and Francis, Philadelphia, 1999).
[57] B. L Yakobson, C. J. Brabec, and J. Bernholc, Phys. Rev.
Lett. 76, 2511 (1996).